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SIGNED DISTANCE FUNCTION AND ITS DERIVATIVES 
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Abstract. We present a fast convolution-based technique for computing an approximate, signed 
Euclidean distance function at a set of 2D and 3D grid locations. Instead of solving the non-linear 
static Hamilton- Jacobi equation (||V»S'|| = 1), our solution stems from solving for a scalar field (j) 
in a linear differential equation and then deriving the solution for S from its exponent. In other 
words, when S and (j) are related by ^ = exp (^~~^ ^ind (j) satisfies a specific linear differential 
equation corresponding to the extremum of a variational problem, we obtain the Euclidean distance 
function S = — rlog((/)) in the limit as r — 0. This is in sharp contrast to solvers such as the 
fast marching and fast sweeping methods which directly solve the Hamilton- Jacobi equation by 
the Godunov upwind discretization scheme. Our linear formulation provides us with a closed-form 
solution to the approximate Euclidean distance function expressed as a discrete convolution, and 
hence efficiently computed by the Fast Fourier Transform (FFT). Moreover, the solution circumvents 
the need for spatial discretization of the derivative operator thereby providing highly accurate results. 
As r — 0, we show the convergence of our results to the true solution and also bound the error 
for a given value of r. The differentiability of our solution allows us to compute — using a set of 
convolutions — the first and second derivatives of the approximate distance function. In order to 
determine the sign of the distance function (defined to be positive inside a closed region and negative 
outside), we compute the winding number in 2D and the topological degree in 3D and again explicitly 
show that these computations can be performed via fast convolutions. We demonstrate the efficacy 
of our method through a set of experimental results. 
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1. Introduction. Euclidean distance functions (more popularly referred to as 
distance transforms) are widely used in image analysis and synthesis [12j. The task 
here is to assign at each grid point a value corresponding to the Euclidean distance 
to its nearest neighbor from a given point-set. Formally stated: given a point-set 
Y = {Yfc G M^, k G {1, . . . , K}} where D is the dimensionality of the point-set and a 
set of equally spaced Cartesian grid points X, the Euclidean distance function problem 
requires us to assign 

S{X) = imn\\X-Yu\\ (1.1) 

k 

where ||.|| represents its Euclidean norm. In computation geometry, this is the Voronoi 
problem [6] and the solution S{X) can be visualized as a set of cones with the centers 
being the point-set locations {Yk}. The Euclidean distance function problem is a 
special case of the eikonal equation where the forcing function is identically equal to 
1 everywhere and hence satisfies the differential equation 

IIV^II = 1, (1.2) 

barring the point-set locations and the Voronoi boundaries where it is not differen- 
tiable. Here = {Sx-, Sy) denotes the gradients of S. This is a nonlinear differential 
equation and an example of a static Hamilton- Jacobi equation. 
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Since the advent of the fast marching method fi^ ^J, the hterature is replete 
with pioneering works which have successfuhy tackled this problem. Fast marching 
is an elegant technique which solves for S in 0{N log N) time at the given N grid 
locations using the Godunov upwind discretization scheme. Faster methods like the 
fast sweeping method [20| employs Gauss-Seidel iterations and solves for S in 0{N). 
The ingenious work in [T9J gives an 0{N) implementation of the fast marching method 
with a cleverly chosen untidy priority queue data structure. Fast sweeping methods 
have also been extended to the more general static Hamilton- Jacobi equation [11] and 
also for the eikonal equation on non-regular grids [9"^T0]. A Hamiltonian approach to 
solve the eikonal equation can be found in [16j. 

The intriguing aspect of our method is that a nonlinear Hamilton-Jacobi equation 
is obtained in the limit as r ^ of a linear equation. When the distance function S is 
expressed as the exponent of a scalar field specifically (f){X) = exp ^ ~'^^^) ^ ^ and if 

is the solution to a variational problem satisfying its corresponding linear Euler- 
Lagrange equation, we show that as r ^ 0, 6* satisfies Equation 11.21 Consequently, 
instead of solving the non-linear Hamilton-Jacobi equation, one can solve for the func- 
tion (j) (taking advantage of its linearity), and then compute an approximate S for 
small values of r. This computational procedure would be approximately equivalent 
to solving the original Hamilton-Jacobi equation. Our linear approach results in a 
closed-form solution, which can be expressed as a discrete convolution and computed 
in O(A^logA^) time using a Fast Fourier Transform (FFT) [3j where N is the number 
of grid points. This is major advantage of our method as the closed- form solution cir- 
cumvents the need for spatial discretization of the derivative operator in Equation [L2l 
a problem that permeates the Hamilton-Jacobi solvers [131 [T21 [20] • This accounts for 
improved accuracy of our technique. However, a minor caveat of our method being 
that the resultant Euclidean distance function is an approximation since it is obtained 
for a small but non-zero value of r, but nevertheless converges to the true solution as 
r ^ 0. 

The linear approach gives us an unsigned distance function. We complement 
this by independently finding the sign of the distance function in 0(7VlogA^) time 
on a regular grid in 2D and 3D. We achieve this by computing the winding number 
for each location in the 2D grid and its equivalent concept, the topological degree in 
3D. We show that just as in the case of the unsigned Euclidean distance function, the 
winding number and the topological degree computations can also be written in closed- 
form, expressed as discrete convolutions and efficiently computed using FFTs. We are 
not cognizant of any previous work that uses the winding number and topological 
degree approaches to compute signed distance functions. Very often, we also seek 
the gradient of the signed distance function which are not easy to obtain via the 
Hamilton-Jacobi solvers due to the lack of differentiability of their solution. Since our 
method results in a differentiable closed-form solution, we can leverage it to compute 
these quantities. Surprisingly we noticed that the solution for the gradients too can be 
written as discrete convolutions and efficiently computed using FFTs. Furthermore, 
they can also be shown to converge to their true value as r ^ 0. To our knowledge, 
fast computation of the derivatives of the distance function on a regular grid using 
discrete convolutions is new. 

The paper is organized as follows. In Section [2] we derive the linear equation 
formalism for the Euclidean distance function problem and show convergence of our 
solution to the true solution as r ^ . We provide an approximate closed-form solu- 
tion to compute the distance function and give an error bound between the computed 
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and true distance functions for a given value of r. Section [3] demonstrates how the 
closed-form solution can be represented as a convolution and thus efficiently computed 
using fast Fourier transforms. In Sections H] and El we compute the winding number 
(in 2D), the topological degree (in 3D) and the derivatives of the distance function, by 
again expressing these quantities using discrete convolutions. In Section [6l we show 
anecdotal evidences for the usefulness of our method by furnishing both experimental 
results and comparisons to standard techniques. We finally conclude in Section [71 

2. Linear equation approach for Euclidean distance functions. Consider 
a variational problem 

m) =r' I Wm'dX + j\4>- 4>l?dX (2.1) 

where (p^ — a function of r — represents the intial wave front concentrated around the 
source locations {Yk}^^i in the as r ^ 0. We define (t)J){X) as 

K 

4>l{X) = Y,4>l{X) (2.2) 

k=l 

where (l)\{X) is chosen such that it is square integrable to 1 with its support sequen- 
tially converging towards the point source as r approaches zero and asymptot- 
ically behaves like the square-root of a J function centered around Y^. The square 
integrability to 1 constraint changes its functional form in accordance with the spatial 
dimension, as explicated in the subsequent sections. 

The Euler-Lagrange equation corresponding to the extremum of computed 
over the scalar field (j) is given by the linear equation 

-T^\7^ + ^ = ro, (2.3) 

where V stands for the Laplacian operator. We may be tempted to replace in 
Equation 12.31 with a combination of delta functions each centered around Yk and 
obtain an inhomogeneous screened Poisson partial differential equation. But since 
the delta functions are not square-integrable, they cannot be incorporated into the 
variational framework given in Equation 12. 1[ Defining as in Equation 12.21 forces it 
to behave like the square-root of a. S function as r ^ and hence is square-integrable 
for all values of r. 

To our amazement, we realized that when we relate (/>(X) = exp (^ ~^^^^ ^ and 
(/) satisfies Equation 12. 3[ S asymptotically satisfies the Hamilton- Jacobi equation 11.21 
in the limit r ^ 0, as elucidated under Section 12. 2[ This relationship motivates us 
to solve the linear equation 12.31 instead of the non-linear eikonal equation and then 
compute the distance function via 



S{X) = -rlog^{X). (2.4) 

2.1. Solution for the Euclidean distance function. We now derive the so- 
lution for (/)(X) (in ID^ 2D and 3D) satisfying Equation 12.31 and then for S{X) from 
the relation [2. 4[ 

Since it is meaningful to assume that S{X) approaches infinity for points at 
infinity, we can use Dirichlet boundary conditions ^(X) = at the boundary of an 
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unbounded domain. Using a Green's function approach [2j, we can write expressions 
for the solution (j). The Green's function G satisfies the relation 

(-r2v' + l)G(X) = -(5(X). (2.5) 

The form of the solution for G [2j in ID, 2D and 31) over an unbounded domain with 

vanishing boundary conditions at oo is given by, 

ID: 

G(X) = ^exp(^^). (2.6) 

2D: 

GiX) = ^,Kjm (2.7) 



27rr2 ^ 
-K"^) ll^l 



2rv/2^r||X|r 



>0.25 



where Kq is the modified Bessel function of the second kind. 
3D: 

The solution for (j) can then be obtained via convolution as 

K 

<P{X) = G{X) * 4>UX) = J2 G(^) * <^fe(^) (2-9) 

k=l 

from which S can be recovered using the mathematical relation 12.41 

2.2. Proofs of convergence. We now prove that as r ^ 0, S{X) — obtained 
from the exponent of (j) which satisfies Euler-Lagrange equation refeqiphiEquationwithtau- 
converges to the true solution r{X) = min/c ||X — Yk\\ = \\X — Yk^W. We show this 
explicitly for each spatial dimension. 

ID: Since we require (/>^(X) to behave like a square-root of the 5 function centered 
at Yfc as r ^ 0, we define it as 



_ / ^ for n - I < X < n + f ; 

otherwise 



Plugging it in Equation 12.91 and using the expression for the ID Green's function, we 
solve for ^ as 

^(^) = E / exp ^ Uz. (2.10) 

Using the relation (|2.4p . the Euclidean distance function is given by 

S{X) = Or - r log f V / exp ( ^] dz] , (2.11) 



where = rlog(2) + |rlog(r) and the integration region BJ^{X) equals 



Bl{X)= \x-Yk-^-,X-Yk+^- 

In order to show convergence to r{X) as r approaches zero, we define 

_ r 1 if. 
"'= = 1-1 if. 

for each k. Then for sufficiently small t we get 
S{X) < -TlogjTexp 

= a-Tlog(T) 



(2.12) 



X>Yk; 
X<Yu 



\X - n„ + ak, 



(2.13) 



as \X-Y^,^OL^,\\ > \Z\ for ZeBl^. 

On the other hand, since \X — Yj^ — a^^l < \Z\^\/Z G B^ at small values of r, we 
also get 



K 



S{X) >a-rlog<r^exp 



k=l 



\X-Yk-ak^\ 



> Cr — T log(r) — r log < K exp 



\X - Yfco - a/eot 



^)} 



a-Tl0g(T)-Tl0g(/f) 



X - Yfco - - 



(2.14) 



In order to see why the second step in the above relation holds, consider the two 
scenarios (i) X lies on the Voronoi boundary and (ii) X is not a point on the Voronoi 
boundary. If X doesn't lie on the Voronoi boundary, then exist a neighborhood A/p(X) 
around X such that Vr G Np{X), \Y -Yk^ < \Y -Y^l^k. Since \X-ak^^\ e Nplx) 
for sufficiently small values of r, the aforementioned relation is true. On the flip side 
if X is a point on the Voronoi boundary, the closest source point Yk is not uniquely 
defined. However we can unambiguously choose a source point Y^^ and a tq such that 
forr G (0,ro], |X - a^.r - F^J < |X - a^r - F^l, Vfc and |X - F^J < {X-Ykl^k. 
This observation ascertains the above inequality. 

Since Cr, rlogr and rlogi^^ approach zero as r 
\X-Yk,\=riX). 



0, we have limr^o S{X) = 



2D: Let the grid location X and the point source Yk be represented by (x^y) and 
(xk^Vk) respectively. We define (/>^(X) as 

{~ 1 ^k 2" — — ^k H~ "2 1 
yk-^<y<yk^i; 
otherwise 

such that it is square integrable to 1 and behaves like the square-root of the 5 function 
centered at Yk as r tends to zero. Using the expression for the 2D Green's function 
and the relation \2A\ the Euclidean distance function is given by 

Six) = Cr - Tlog (j2^ 1^^^^^ Ko (^^) dZ^ , (2.15) 
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where = Tlog(27r) + 3r logr and the integral region BJ^{X) equals 

r T T 1 r 1 



T r ' 



V-Vk-^,V-Vk^^ 



(2.16) 



Defining 



and 



I -1 i 



if X > Xjt; 
if X < 



if 2/ > Vk\ 
if 2/ < 



for each k and closely following the arguments illustrated for the ID case, we get 
5(X)<a-rlog{.i.o(^i^)} 

= a - r log(r) - rlog ji^o (HlL^) 



(2.17) 



for small values of r, where Xi^- = (x + ce/c^, 2/ + /^fci)- As in the ID case, note that 
ll^ix-lfcoll > ll^ll for^eS,„.^ 

Using the relation Kq,{z) > when z > 0.5, we observe that 



5(X)<a-rlog(r)+rlog. 



\\Xir-Yk,\ 



\Xlr - ifeo 



(2.18) 



as T approaches zero. 

If we let = — <^/c§5 ^ — /^/cf ) 5 then similarly to the ID case we arrive at 
the inequality 



S{X)>Cr-TlogiTj2Ko ( 
I fe=l ^ 



\\X2r-Yk\ 



>Cr-TlogiT)-T\og{KKo 



n\X2r-Yk,\ 



= Cr - rlog(r) - rlog(Jr) - rlog < Kq 



\Xo^- 



^)} 



As Kq{z) < exp(— when 2: > 1.5, we get 

S{X) > Cr - rlog(r) - r log(K) + \\X2r - Y^, 



(2.19) 



(2.20) 



at small values of r. Since Xir,X2T approach X as r ^ and the rest of the terms 
tend to zero in the limit, we arrive at our desired result, namely limT-^o = 
\\X-YkJ=r{X). 

3D: Let the grid location X and the source point Yk be denoted by X = (x^y^z) 
and Yk = {xk^Vk^^k) respectively. To ensure square-integrability to 1 we define 0^ as 



/ 1 . 



f>l{x) 



5 Xk n ^ X ^ Xk H~ 9 , 

Vk - ^ <y <yk-\-^, 

Zk - ^ < z < Zk -\- 
otherwise 



By exactly following the line of argument described for the ID and the 2D case 
where we bound S{X) above and below by functions which converge to the true 
Euclidean distance function r{X) as r tends to zero, we can establish the result 
limr^oS{X) = \\X-Yk,\\=r{X). 

2.3. Closed- form solution and the error bound between the obtained 
and true Euclidean distance function. Based on the nature of the expression 
for the Green's function, it is worth accentuating the following very important point. 
Observe from above that the expression for Green's function G in either ID^ 2D or 
3D takes the form 



exp 



lim ' = 0, for ||X|| ^ (2.21) 

r^o crd\\X\\P ' N N ^ V / 

for c, d and p being constants greater than zero. In the limiting case of r — ^ 0, we see 
that if we define 

G{X) = Cexp (^^^) , (2.22) 

for some constant C, then 

lim \G{X) - G{X)\ = 0, for ||X|| ^ (2.23) 

and furthermore, the convergence is uniform for ||X|| away from zero. Therefore, 
G{X) provides a very good approximation for the actual unbounded domain Green's 
function as r ^ 0. For a fixed value of r and X, the difference between the Green's 

functions is O ( — — - ) which is relatively insignificant for small values of r and 

for all X 7^ 0. Moreover, using G also also avoids the singularity at the origin in the 
2D and 3D case. The above observation motivates us to compute the solutions for (j) 
by convolving with G instead of the actual Green's function G. 

Furthermore, our proof technique used to manifest convergence of 5'(X) to the 
true solution r(X) — described in the preceding section — also encourages us to sup- 
plant the integral f]gT(^x) G{Z)dZ — obtained as a result of convolving the Green's func- 
tion G with 0^ — with r^G{X — Yj.) at small values of r. Here D corresponds to the 
spatial dimension. This is a cogent approximation for the following reasons. Firstly, 
the integral region S^(X) is considered around X — Y^ with its area/volume dwindling 
to zero as r approaches zero (refer Equations 12.121 and I2.16p . Hence we can replace 
the integral JjgTf^x) G{^)dZ with its Reimann summation by assuming that G{Z) is 
constant over S^(X) and equals the mid-point value G{X — Yfc). Secondly, our proof 
explicitly corroborates that substituting the integral by either sup^^23^(x) ^(^) 
or infzG23^(x) G{Z) , we can still establish convergence to the true solution. Since 
i'nfzeBi{x)G{Z) < G{X — Y^) < ^^VzeB^{x) G{Z)^ our approximation is sound at 
small values of r. Thirdly, it bestows upon us with a closed-form solution for (j) as 
seen below. 

These insights inspires us to solve for the scalar field (j) as 

HX) « |r^r--C g exp (^^;^) | (2.24) 
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rather than via Equation 12.91 Here corresponds to the value of (j)^{Yk)^ with 7 
determined by the spatial dimension. The approximate Euclidean distance function 
computed based on the relation 12.41 is given by 

S{X) « (7 - D)t log t - t log C - t log 1^ exp | . (2.25) 

Since (7 — D)r log r and r log C are constants independent of X and converges to zero 
as r ^ 0, they can be ignored while solving for S at small values of r. Hence the scalar 
field ^(X) corresponding to the Euclidean distance function can be approximated by 

</.(X)^f:exp(^^^). (2.26) 

and the approximate Euclidean distance function equals 

S{X) = -r log I f exp (^d^^) | . (2.27) 

We would like to underscore that the approximate function defined in Equation 12.261 
possess all the desirable properties that we need. Firstly, we notice that as r ^ 
0, (l){Yk) ^ 1 at the given point-set locations Y^. Using the relation ()2.4p we get 
S{Yk) ^ as r ^ satisfying the initial conditions. Secondly for small values of 

r, ^^P ^ can be subrogated by exp where r{X) is the true 

Euclidean distance given by r{X) = min/^ \\X — 1/^11- Hence 



S{X) ^ -r log exp ^ ^ J = r(X). (2.28) 

Thirdly, (/>(X) can be easily computed using the fast Fourier transform as discussed 
under the subsequent section. Hence for all practical, computational purposes we 
consider the function shown in Equation 12.261 The bound derived below between the 
approximate S{X) and the true solution r{X) also unveils the proximity between the 
computed and the actual Euclidean distance function. 
Note from Equation 12.271 that 

S{X) < -r log exp (^^^^) = r{X). (2.29) 

Also we get 

S{X) > -r log ji^exp (^^^^) } 

= -r log i^ + r(X) (2.30) 

and hence, 

r{X)-S{X) <r\ogK. (2.31) 

From Equations 12.291 and 12. 3H we have 

\r{X)-S{X)\ <T\ogK. (2.32) 

It is worth commenting that the bound r log K is actually very tight as (i) it scales 
only as the logarithm of the cardinality of the point-set (K) and (ii) it can be made 
arbitrarily small by choosing a small but non-zero value of r. 



Table 3.1 

Approximate Euclidean distance function algorithm 

~1. Compute the function f{X) = exp ^^^^^ at the grid locations. 

2. Define the function Skron{X) which takes the value 1 at the point-set locations 
and at other grid locations. 

3. Compute the FFT of / and namely Ffft{U) and Gfft{U) respectively. 

4. Compute the function H{U) = Ffft{U)Gfft{U). 

5. Compute the inverse FFT of H{U) to obtain ^(X). 

6. Take the logarithm of ^(X) and multiply it by (— r) to recover 
the approximate Euclidean distance function. 



3. Efficient computation of the approximate unsigned Eucfidean dis- 
tance function. We now provide a fast 0{N log N) convolution-based method to 
compute the distance transform on a set of N grid locations {Xi^i = {1, . . . ^N}}. 
The solution for (j){X) in Equation 12.261 at the grid locations can be represented as 
the discrete convolution between the functions 

/(X) ^ exp ( ^) (3.1) 



computed at the grid locations, with the function g{X) which takes the value 1 at the 
point-set locations and at other grid locations, i.e, 

K 

g{X) = Y,hron{X-Yk) (3.2) 

k=l 

where, 

^^roniX-n)^[l 11=^^^ (3.3) 

By the convolution theorem [3j, a discrete convolution can be obtained as the 
inverse Fourier transform of the product of two individual transforms, which for two 
0{N) sequences can be computed in 0{N log N) time [5j. One just needs to compute 
the discrete Fourier transform (DFT) of the sampled version of the functions f{X) 
and g{X), compute their point- wise product and then compute the inverse discrete 
Fourier transform. Taking the logarithm of the inverse discrete Fourier transform 
and multiplying it by (— r), gives the approximate Euclidean distance function. The 
algorithm is adumbrated in Table 13.11 

3.1. Computation of the approximate Eucfidean distance function in 
higher dimensions. Our technique has a straightforward generalization to higher 
dimensions. Regardless of the spatial dimension, the approximate Euclidean distance 
function, S can be computed by exactly following the steps delineated in Table 13. li 
It is worthwhile mentioning that computing the discrete Fourier transform using the 
FFT is always 0{N logN) irrespective of the spatial dimension. Hence, for all dimen- 
sions, S can be computed at the given N grid points in 0(7Vlog7V). This speaks for 
the scalability of our technique, which is generally not the case with other methods, 
for example KD-Trees [6j B. 



* Though the actual number of grid points (A/") increases with dimension, the solution is always 
0(N log N) in the number of grid points. 
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3.2. Numerical issues. In principle, we should be able to apply our technique 
at very small values of r and obtain highly accurate results. But we noticed that a 
naive double precision-based implementation tends to deteriorate for r values very 
close to zero. This is due to the fact that at small values of r, f{X) drops off 
very quickly and hence for grid locations which are far away from the point-set, the 
convolution done using FFT may not be accurate. To this end, we turned to the 
GNU MPFR multiple-precision arithmetic library which provides arbitrary precision 
arithmetic with correct rounding [7J. MPFR is based on the GNU multiple-precision 
library (GMP) [T8]. It enabled us to run our technique at very small values of r giving 
highly accurate results. We corroborate our claim and demonstrate the usefulness of 
our method with the set of experiments described under Section [6l 

3.3. Exact computational complexity. More the number of precision bits p 
used in the GNU MPFR library, better is the accuracy of our technique, as the error 
incurred in the floating point operations can be bounded by 0(2~^). But using more 
bits has an adverse effect of slowing down the running time. The 0{N log N) time 
complexity of the FFT algorithm [5j for an 0{N) length sequence represents only the 
number of floating-point operations involved, barring any numerical accuracy. The 
accuracy of the FFT algorithm and our technique entirely depends on the number of 
precision bits used for computing elementary functions like exp, log, sin and cos and 
hence should be taken into account while calculating the exact time complexity. If p 
precision bits are used, these elementary functions can be computed in 0{M{p) logp) 
[31|lll[T7], where M{p) is the computational complexity for multiplying two p-digit 
numbers. The Schonhage-Strassen algorithm [15j gives an asymptotic upper bound 
for M{p) with a run-time bit complexity of M{p) = 0{p logp log logp). The actual 
running time of our algorithm — while taking these p precision bits into account — 
for computing S at the given N grid locations is then 0{N log{N)p{logp)^ log(logp)) 
bit- wise operations. 

4. Fast computation of signed distance functions. The solution for the 
approximate Euclidean distance function in Equation 12.271 is lacking in one respect: 
there is no information on the sign of the distance. This is to be expected since the 
distance function was obtained only from a set of points Y and not a curve or surface. 
We now describe a new method for computing the signed distance function in 2D 
using winding numbers and in 3D using the topological degree. 

4.1. Computing winding numbers. Assume that we have a closed, paramet- 
ric curve {x^^\t) ^ x^'^\t)} , t G [0, 1]. We seek to determine if a grid location in the 
set {Xi gM^, i G {!,..., N}} is inside the closed curve. The winding number is the 
number of times the curve winds around the point Xi (if at all) with counterclockwise 
turns counted as positive and clockwise turns as negative. If a point is inside the 
curve, the winding number is a non-zero integer. If the point is outside the curve, 
the winding number is zero. If we can efficiently compute the winding number for 
all points on a grid w.r.t. to a curve, then we would have the sign information (in- 
side/outside) for all the points. We now describe a fast algorithm to achieve this 
goal. 

If the curve is C^, then the angle 6{t) of the curve is continuous and differentiable 
and dO{t) = y^p^ ^ ^ dt. Since we need to determine whether the curve winds 

around each of the points Xi,i G {!,..., N}, define x-^^) = (x*^-^^ — x*^^^ — 
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Vz. Then the winding numbers for the grid point Xi is 



dt, \/ie{l,...,N}. (4.1) 



As it stands, we cannot actuahy compute the winding numbers without performing 
the integral in Equation 14. 1[ To this end, we discretize the curve and produce a 
sequence of points {Y^ G M^, A: G {1, . . . , K}} with the understanding that the curve 
is closed and therefore the "next" point after Yk is Yi . (The winding number property 
holds for piecewise continuous curves as well.) The integral in Equation 14. II becomes 
a discrete summation and we get 



Vz G {1, . . . , A^}, where the notation Y^^i denotes that Y^'^-^ = Y^'_^-^ for G {1, . . . , K— 

1} and Y^0i =Y^'\ We can simplify the notation in Equation l4.2l fand obtain a mea- 
sure of conceptual clarity as well) by defining the "tangent" vector {Z^^k = {1,..., K}} 

as 

zi^=Y^2i-Yl-\k€{l,...,K} (4.3) 

with the (•) symbol indicating either coordinate. Using the tangent vector Z/e, we 
rewrite Equation 14.21 as 

I « (n'" -^.'") - (if - xP) 4" , , , , 

We now make the somewhat surprising observation that /i^ in Equation 14.41 is 
a sum of two discrete convolutions. The first convolution is between two functions 
fcr{X) = fc{X)fr{X) and g2{X) = Ylk=i ^^^''^kroniX -Yk). The second convolution 
is between two functions fsr{X) = fs{X)fr{X) and ^i(X) = Xl^Li zj^'' ^kron{X -Yk). 
The Kronecker delta function 6kron{^ — Yk) is defined Equation 13.31 The functions 
fc{X), fs{X) and fr{X) are defined as 

/,(X) ^ MX) ^ and (4.5) 

MX) . ^ (4.6) 

with the understanding that /c(0) = /s(0) = /r(0) = 0. Here we have abused notation 
somewhat and let X*^^^ (X*^^^) denote the x (?/)-coordinate of the grid point X. Armed 
with these relationships, we rewrite (j4.4p to get 

KX) = ^ [-fcr{X)*g2{X)+U{X)*g,{X)] (4.7) 

which can be computed in O(A^logA^) time using FFT-based convolution simultane- 
ously for all the N grid points {X^, i = {1, . . . , X}}. 
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4.2. Computing topological degree. The winding number concept for 2D 
admits a straight forward generahzation to 3D and higher dimensions. The equivalent 
concept is the topological degree which is based on normalized flux computations. 
Assume that we have an oriented surface in 3D pj which is represented as a set of 
K triangles. The k^^ triangle has an outward pointing normal Pk and this can easily 
be obtained once the surface is oriented. (We vectorize the edge of each triangle. 
Since triangles share edges, if the surface can be oriented, then there's a consistent 
way of lending direction to each triangle edge. The triangle normal is merely the 
cross-product of the triangle vector edges.) We pick a convenient triangle center (the 
triangle incenter for instance) for each triangle and call it Y^. The normalized flux 
(which is very closely related to the topological degree) [I| determines the ratio of the 
outward flux from a point Xi treated as the origin. If Xi is outside the enclosed surface, 
then the total outward flux is zero. If the point is inside, the outward normalized flux 
will be non-zero and positive. 

The normalized flux for a point Xi is 

This can be written in the form of convolutions. To see this, we write Equation 14.81 
in component form: 

which can be simplified as 

/^W = (/i(^) * ^i(^) + h{X) * g2{X) + h{X) * g^{X)) (4.10) 

where f(.){X) = and g(.){X) = ^fe'^4ron(^ - Yk) where the Kronecker 

delta function Skroni^ — ^k) is defined Equation 13.31 This can be computed in 
0{N log N) time using FFT-based convolution for all the N grid points X^. 

For the sake of clarity we explicitly show the generalization of the winding number 
to the topological degree by rewriting some of the calculations involved in computing 
the winding number. Recall that for every point Yk on the discretized curve, we 
defined its tangent vector Zk as in equation (j4.3p . The outward pointing normal 
Pk = (P^^\ P^^^), at the point Yk {Pk will point outwards provided Yi, I2, • • • , Y/e are 
taken in the anti-clockwise order), is given by P^^^ = Z^\p^^ = —Z^\ Using the 
normal vector P/e, equation (|4.4p can be rewritten as 

u. = ^f«^^^^^ (411) 

Notice the similarity between equations (j4.1ip and (j4.8p . This manifests that the 
topological degree is just a generalization of the winding number concept. 

We have thus demonstrated that the sign component of the Euclidean distance 
function can be separately computed (without knowledge of the distance) in parallel 
in 0(A^log7V) on a regular 2D and 3D grid. 
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5. Fast computation of the derivatives of the distance function. Just 
as the approximate Euclidean distance function S{X) can be efficiently computed in 
0{N log A/"), so can the derivatives. This is important because fast computation of the 
derivatives of S{X) on a regular grid can be very useful in curvature computations. 
Below, we detail how this can be achieved. We begin with the gradients and for 
illustration purposes, the derivations are performed in 2D: 

SAX) = ; ), (5.1) 

Ef=.exp(-M^) 

A similar expression can be obtained for Sy{X). These first derivatives can be rewrit- 
ten as discrete convolutions: 



^^^""^ - f{X)*9{X) ' - f{X).g{X) ' 



where fc{X) and fs{X) are as defined in Equation 14.51 and f{X) and g{X) are given 
in Equations 13.11 and 13.21 respectively. 

The second derivative formulas are somewhat involved. Rather than hammer out 
the algebra in a turgid manner, we merely present the final expressions — all discrete 
convolutions — for the three second derivatives in 2D: 



T 



Q fY^ [-7fsiX) + fAX)fr{X)]f{X)*g{X) , 1,,,,,^, , 
Syy{X) = — ^ f{X)*g{X) ^ ^' ^ 



[7 + fAX )]f,{X)f,iX)fiX)*g{X) 1 
f{X)*g{X) 



SccyiX) = 'Y/', !U + -SAX)Sy{X) (5.5) 



where fr{X) is as defined in Equation! 

Since we can efficiently compute the first and second derivatives of the approxi- 
mate Euclidean distance function everywhere on a regular grid, we can also compute 
derived quantities such as curvature (Gaussian, mean and principal curvatures) for 
the two-dimensional surface S{X) computed at the grid locations X. In Section [6l 
we visualize the derivatives for certain shape silhouettes. 

5.1. Convergence analysis for the derivatives. We now show convergence 
of the derivatives (Sx^Sy) obtained via Equation 15.11 to their true value as r ^ 0. 
Recall that the true distance function is not differentiable at the point-source locations 
{Yk}^^i and on the Voronoi boundaries which corresponds to grid locations which 
are equidistant from two or more point sources. Hence it is meaningful to speak 
about convergence only for grid locations whose closest source point can be uniquely 
determined. 

For illustration purposes we show the analysis in 2D. Let Yk^ denote the unique 
closest source point for the grid location Xq, i.e, \\YkQ — -^o|| < \\Yk — Xo||, V/c / ko. 
Then its true derivatives are given by 

^^^""^^ ^ Ilio-i^oV^^*^''^^ ^ \\Xo - • ^'-'^ 
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Multiplying and dividing the Equation 15. II by exp 

v(^) \-(^) ^(1) ,.(1) 

11-^0-^0 II +^fc=l,Mfco ll^o-nil ^^Py r) 
^x[Xo) = —- (5.7) 

where = ||-^o — ^fcll — ||-^o — ^/eoll- Since > 0,V/c 7^ /cq, it follows that 
lim^-^o exp (— ^) = 0,V/c 7^ /cq. Since all the other terms in Equation 15.71 are inde- 
pendent of r, we get liniT-^o ^'^^(-^o) = 5'*(Xo). The convergence analysis for Sy{Xo) 
follows along similar lines. 

Furthermore, the gradient magnitude ||VS'(X)|| for any non-zero value of r will 



be strictly less than its true value y + = 1. To see this consider 

Equation 15.11 and for the sake of convenience define = "pr^VklT^ ~ ||j^_y^^|| 
and = exp (^-^^^^^ Since a| + = 1, Vfc, we get 

llV7cr^M|2 Ef^i -^1 + 2 Ef^i E^fe (o^fco^^ + &k&i) EkEi 

l^k=l ^fe + ^ l^k=l l^l>k ^kEl 

Using the Cauch-Schwarz inequality, we have 

WkO^i + Mi\ < ^o^l^Pl^o^.^pf = 1. (5.9) 

It is then easy to see that ||V6'(X)|| < 1. We provide experimental evidence in 
the subsequent section to corroborate this fact. Nevertheless, the magnitude value 
converges to 1 as r ^ at all the grid locations barring the point-sources and the 
Voronoi boundaries, as the gradients themselves converge to their true value. 

6. Experiments. In this section we show the usefulness of our linear approach 
to computing Euclidean distance functions by running it on a bounded 2D and 3D 
grid. As we discussed before, in order improve the computational accuracy of our 
technique we are forced to go beyond the precision supported by the double fioating- 
point numbers (64 bits) and resort to arbitrary precision packages like GNU multiple- 
precision library (GMP) [18j and MPFR [7J. For the following experiments we used 
p = 512 precision bits. 

6.1. 2D Experiments. Example l:We begin by demonstrating the effect of 
r on our method and evince that as r ^ 0, the accuracy our method does improve 
significantly. To this end, we considered a 2D grid consisting of points between 
(-0.121,-0.121) and (0.121,0.121) with a grid width of ^. The total number of 
grid points is then TV = 125 x 125 = 15,625. We ran 1000 experiments each time 
randomly choosing 5000 grid locations as data points (point-set), for 9 different values 
of r ranging from 5 x 10~^ to 4.5 x 10~^ in steps of 5 x 10~^. For each run and each 
value of r, we calculated the percentage error as 

error = ^2^^, (S-l) 

i=l * 

where Di and are respectively the actual distance and the absolute difference 
of the computed distance to the actual distance at the i^^ grid point. Figure 16.11 
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shows the mean percentage error at each value of r. The maximum value of the 
error at each value of r is summarized in table 16.11 The error is less than 0.6% at 
r = 0.00005 demonstrating the algorithm's ability to compute accurate Euclidean 
distance functions. 




% 0.5 1 1.5 2 2.5 3 3.5 4 4.5 
Tau Values .,„-4 



Figure 6.1. Percentage error versus r in 1000 2D experiments. 



r 


Maximum error 


0.00005 


0.5728% 


0.0001 


1.1482% 


0.00015 


1.7461% 


0.0002 


2.4046% 


0.00025 


3.1550% 


0.0003 


4.0146% 


0.00035 


4.9959% 


0.0004 


6.1033% 


0.00045 


7.3380% 



Table 6.1 

Maximum percentage error for different values of r in 1000 2D experiments. 



Example 2: We then executed our algorithm on a set of 2D shape silhouettes [16|1|. 
The grid size is -0.125 < x < 0.125 and -0.125 <y< 0.125 with a grid width of 
2^. The number of grid locations equals N = 257 x 257 = 66,049. We set r for 
our method at 0.0003. For the sake of comparison, we ran the fast sweeping method 
for 10 iterations sufficient enough for it to converge. The percentage error calculated 
according to Equation l6.1l for both our method and fast sweeping in comparison to the 
true Euclidean distance function for each of these shapes are adumbrated in Table [Ql 

The true Euclidean distance function contour plots and those obtained from our 
method and fast sweeping are delineated in Figure [621 

In order to differentiate between the grid locations that are either inside or outside 
each shape, we computed the winding number for all the grid points simultaneously in 
O(A^logA^) using our convolution-based winding number method. Grid points with 
a winding number value greater than after rounding were marked as interior points. 



^We thank Kaleem Siddiqi for providing us the set of 2D shape silhouettes used in this paper. 
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Shape 


Our linear method 


Fast sweeping 


Hand 


2.182% 


2.572% 


Horse 


2.597% 


2.549% 


Bird 


2.116% 


2.347% 



Table 6.2 

Percentage error for the Euclidean distance function computed using the grid points of these 
silhouettes as data points 




Figure 6.2. Contour plots: (i) Left: True Euclidean distance function, (ii) Center: Our 
method, (Hi) Right: Fast sweeping 

In the the left part of Figure [Ql we visualize the vector fields (Sx^Sy) for all the 
interior points (marked in blue). We see that our convolution-based technique for 
computing the winding number separates the interior grid points from the exterior 
with almost zero error. In the right part of Figure 16.31 we plot the distribution of 
the winding number values computed over all the interior and the exterior locations. 
Observe that for almost all the grid points, the winding number values are close to 
binary, i.e either or 1, with a value of marking the exterior points (as the curve 
doesn't wind around them) and 1 representing interior locations. A zoomed version 
of the quiver plot is shown in Figure 16. 4[ 

In the left part of Figure 16. 5[ we picturize the distribution of the gradient magni- 
tude (IIV^II) values. Since we don't solve for the true distance function S and rather 
approximate it by Equation I2.27[ the magnitude of the gradients (Sx^Sy)) obtained 
via Equation 15. II mav not necessarily satisfy the Hamilton- Jacobi equation II. 2[ Their 
distribution values portrayed in Figure 16.51 manifests it. The reformulation of the 
derivatives given in Eg nation 15 . 71 drives this point home. Nevertheless, we do observe 
that most of the gradient magnitude values (about 90%) are greater than 0.9. Also 
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Distribution of winding number vaiues 






Winding number 
Distribution of winding number values 



Winding number 
Distribution of winding number values 



Winding number 
Distribution of winding number values 



Winding number 



Figure 6.3. (i) Left: Quiver plot ofWS = {Sx,Sy) for a set of silhouette shapes (best viewed 
in color), (ii) Right: Distribution of winding numbers 



notice that all the gradient magnitudes are less than 1, as we evinced in Section ISTTl 
In the right part of Figure 16. 5[ we visualize the gradient magnitude values as an 
image plot in tints of gray. Grid locations whose gradient magnitude values are in 
proximity to 1 are marked as white and black indicates grid points whose gradient 
magnitude values are closer to its minimum value on the grid. From these image plots 
it is quite clear that the grid points which are either too close to the point-sources or 
lie along the medial axis (corresponding to the Voronoi boundaries for these shapes) 
incur maximum errors in their gradient magnitude values. As these grid locations are 
almost equidistant from multiple point-sources, many of them contribute substantially 
to the summation value exp ^_ ll-^-^fcll ^ instead of just the nearest one and hence 
the induced error is relatively high. Interestingly, it might be possible to actually 
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Figure 6.4. Zoomed quiver plot 



leverage these inaccuracies for something good as they might aid in the medial axis 
computation. Notice that, it is fairly an easy task to trace the medial axis from these 
image plots. We are currently investigating whether this can be achieved. 

6.2. 3D Experiments. Example 3:We also compared our Euclidean distance 
function algorithm with the fast sweeping method [2Qj and the exact Euclidean dis- 
tance on the ^^Dragon^^ point-set obtained from the Stanford 3D Scanning Repos- 
itory!. The common grid was -0.117 < x < 0.117, -0.086 < y < 0.086 and 
-0.047 < z < 0.047 with a grid width of ^. We ran our approach at r = 0.0004 
and ran the fast sweeping method for 15 iterations which is sufficient for the Gauss- 
Seidel iterations to converge. We then calculated the percentage error for both these 
techniques according to Equation 16.11 While the average percentage error for our 
approach when compared to the true distance function was just 1.306%, the average 
percentage error in the fast sweeping method was 6.84%. Our FFT-based approach 
does not begin by discretizing the spatial differential operator as is the case with the 
fast marching and fast sweeping methods and this could help account for the increased 
accuracy. 

The isosurface obtained by connecting the grid points at a distance of 0.005 from 
the point set, determined by the true Euclidean distance function, our algorithm and 
fast sweeping are shown in Figure 16.61 The similarity between the plots provides 
anecdotal visual evidence for the usefulness of our approach. 

Example 4: Next, to demonstrate the efficacy of our convolution-based technique for 
computing the topological degree, we ran the following experiments in 3D. The grid 
was confined to the region -0.125 < x < 0.125, -0.125 <y< 0.125 and -0.125 < 
z < 0.125 with a grid width of ^. The number of grid points was N = 274,625. 
Given a set of points sampled from the surface of a 3D object, we triangulated the 
surface using some of the built-in MATLAB® routines. We consider the incenter of 
each triangle to represent the data points {Yk}^^^. The normal Pk for each triangle 



^This dataset is available at |http:/ /graphics. stanford.edu/data/3Dscanrep/[ 
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Figure 6.5. Gradient magnitude (\\VS\\) (i) Left: Distribution, (ii) Right: Image plot 



Figure 6.6. Isosurfaces: (i) Left: Actual Euclidean distance function, (ii) Center: Our algo- 
rithm and (Hi) Right: Fast sweeping 



can be computed from the cross-product of the triangle vector edges. The direction 
of the normal vector was determined by taking the dot product between the position 
vector Yk and the normal vector P^. For negative dot products, Pk was negated to 
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obtain a outward pointing normal vector. We then computed the topological degree 
for all the N grid locations simultaneously in 0{N log N) by running our convolution- 
based algorithm. Grid locations where the topological degree value exceeded 0.7 were 
marked as points lying inside the given 3D object. Figure [6771 shows the interior points 
for the three 31^ objects — cylinder, cube and sphere (left to right). 



Sampled points from the surface Sampled points from the surface Sampled points from the surface 




Figure 6.7. Topological Degree: (i) Top: Sampled points from the surface, (ii) Bottom: Grid 
points lying inside the surface (marked as blue) 



7. Conclusion. In this work, we furnished a linear, variational formalism for 
the Euclidean distance function problem where we computed solutions on a bounded 
2D and 3D grid. We posed a specific variational problem and evinced that the solu- 
tion to its linear Euler-Lagrange equation at small values of r can be used to obtain 
the Euclidean distance function. The intriguing aspect of our approach is that the 
non-linear Hamilton-Jacobi equation is embedded inside a linear equation and the 
solution is derived in the limiting case of r 0. We initially derived the solution 
for (j) satisfying the Euler-Lagrange equation via the Green's function approach and 
later approximated it with a closed- form solution, the major advantage being that the 
closed-form solution is representable as a discrete convolution which can be efficiently 
computed in 0{N\ogN) using FFT. The Euclidean distance is finally recovered from 
its exponent. Since the scalar field (j) is computed for a small but non-zero r, the ob- 
tained Euclidean distance function is an approximation. We derived analytic bounds 
for the error of the approximation for a given value of r and provided proofs of con- 
vergence to the true distance function as r ^ 0. The differentiability of our solution 
endows us to compute the gradients and curvature quantities of the distance function 
S in closed-form, also written as convolutions. Finally, we demonstrated how our dis- 
crete convolution-based technique for computing the winding number in 2D and the 
topological degree in 3D are useful in determining the sign of the distance function. 

While Hamilton-Jacobi solvers have gone beyond the eikonal equation and regular 
grids — by providing efficient solutions even for the more general static Hamilton- 
Jacobi equation on irregular grids [UlElIinj — in the current work we restricts ourselves 
only to computing the Euclidean distance function on regular grids. In the future, we 
would like to follow the pioneering works of the fast marching [13j and fast sweeping 
|2U| methods and try to extend our linear formalism even to irregular grids. 
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